In this project, the open-source R programming language is used to model/monitor changes in the gait cycles over time. R is maintained by an international team of developers who make the language available at The Comprehensive R Archive Network. Readers interested in reusing our code and reproducing our results should have R installed locally on their machines. R can be installed on a number of different operating systems (see Windows, Mac, and Linux for the installation instructions for these systems). We also recommend using the RStudio interface for R. The reader can download RStudio for free by following the instructions at the link. For non-R users, we recommend the Hands-on Programming with R for a brief overview of the software’s functionality. Hereafter, we assume that the reader has an introductory understanding of the R programming language.
In the code chunk below, we load the packages used to support our analysis.
if(require(checkpoint)==FALSE) install.packages("checkpoint") # check to see if checkpoint is installed; if not, install it
library(checkpoint) # package used to facilitate the reproducibility of our work
# a checkpoint of R packages on CRAN on June 08, 2020 to enable the reproduction of our work in the future
checkpoint("2020-06-08")
# check if packages are not installed; if yes, install missing packages
packages = c("tidyverse", "magrittr", # typical data analysis packages
"foreach", "parallel", # packages used for parallelizing some code chunks
"MALDIquant", # match closest points between two vectors
"MASS", # calculates the generalized inverse matrix
"R.matlab")
newPackages = packages[!(packages %in% installed.packages()[,"Package"])]
if(length(newPackages) > 0) install.packages(newPackages)
# using the library command to load all packages; invisible used to avoid printing all packages and dependencies used
invisible(lapply(packages, library, character.only = TRUE))
source("./Functions.R") # our custom built functions
set.seed(2020)
startTime <- Sys.time()
In this paper, we perform a secondary data analysis of the experimental data reported in Baghdadi et al. (2019), where fifteen participants were equipped with an IMU sensor wrapped on their right ankle. We used the IMU data to transform the linear acceleration signals from the local frame of reference to the global frame of reference. Next, we calculated the acceleration magnitude signal and smoothed it by applying a low pass butterworth filter with a cut-off frequency of 10 Hz. Further, we segmented the idividual gait cycles from foot-flat to foot-flat for each participant by tuning the parameters of the segmentation algorithm for each participant.
We load four different sets of data into R:
The code chunk below provides our code for reading the three different datasets. The subject heights dataset is saved into a vector, which we use to scale the computed gait/stride length and height. The main outputs from this chunk are two lists (a) subjectFeaturesList which contains the scaled stride length, scaled stride height, stride duration and experimental time, and (b) subjectProfilesList which contains the sample IMU acceleration points making each acceleration profile.
# Reading the subject Heights from the Baghdadi et al. (2019) GitHub Repo:
# ------------------------------------------------------------------------
subject_heights <- read.csv("https://raw.githubusercontent.com/fmegahed/fatigue-changepoint/master/Data/Raw/DemographicData.csv") %>% pull(Height) %>% `*` (0.01) %>% .[-13] # pull height and convert into meters
# Computing gait features based on experimental data from Baghdadi et al. (2019):
# -------------------------------------------------------------------------------
featureFiles <- list.files(path = "../Data/matFiles/Features/", pattern = "*.mat", full.names = TRUE)
rawSubjectsData <- lapply(featureFiles, readMat) %>% # read all mat files in Feature List
lapply("[[", "M.i.k") # extract sublist titled "M.i.k" from all the lists
names(rawSubjectsData) <- str_extract_all(featureFiles, "Subject[:digit:]*", simplify = TRUE)
numRows <- lengths(rawSubjectsData)/17 # number of observations per Subject
subjectFeaturesList <- vector(mode = "list", length = length(rawSubjectsData)) # initialize list
names(subjectFeaturesList) <- names(rawSubjectsData)
for (i in seq(length(rawSubjectsData))) {
gaitLength <- rawSubjectsData[[i]][ 1:numRows[i] ] %>%
lapply(range) %>% lapply(diff) %>% unlist() %>%
`/` (subject_heights[i])
gaitHeight <- rawSubjectsData[[i]][ (1+numRows[i]):(2*numRows[i]) ] %>%
lapply(range) %>% lapply(diff) %>% unlist() %>%
`/` (subject_heights[i])
gaitDuration <- rawSubjectsData[[i]][ (1+ 2*numRows[i]):(3*numRows[i]) ] %>%
lapply(range) %>% lapply(diff) %>% unlist(recursive = TRUE)
expTime <- rawSubjectsData[[i]][ (1+ 16*numRows[i]):(17*numRows[i]) ] %>%
lapply(unlist) %>% lapply(min) %>% unlist()
subjectFeaturesList[[i]] <- cbind(gaitLength, gaitHeight, gaitDuration, expTime)
}
# Reading the segmented gait cycles for each participant:
# -------------------------------------------------------
profileFiles <- list.files(path = "../Data/matFiles/Profiles/", pattern = "*.mat", full.names = TRUE)
subjectProfilesList <- lapply(profileFiles, readMat) %>% # read all mat files from the path above
lapply("[[", "gait") # extract sublist titled "gait" from all the lists
names(subjectProfilesList) <- sapply(profileFiles, str_extract_all, "Subject[:digit:]*", simplify = TRUE)
# Reading task-related based on the tagging of the experimental videos:
# ---------------------------------------------------------------------
videoFiles <- list.files(path = "../Data/csvFiles/", pattern = "*.csv", full.names = TRUE)
videoDataList <- lapply(videoFiles, read.csv)
names(videoDataList) <- sapply(videoFiles, str_extract_all, "Subject[:digit:]*", simplify = TRUE)
In Section 2, we have read four different data sets and extracted the three gait features. In this subsection, we perform further preprocessing of the data. Specifically, we remove the first 10 minutes of experimental data from the videoDataList since these corresponded to the warm-up period of the experiment as detailed in Maman et al. (2017), Baghdadi et al. (2019), and Maman et al. (2020). WE NEED TO CONTINUE THE DISUCSSION OF WHAT WE HAVE DONE IN THIS SECTION. I STOPPED WORKING ON THIS UNTIL WE DISCUSS. Saeb, this is currently slightly different from your implementation as I am removing any task that spanned the warm-up period; I think this is a cleaner way of doing it and it will not change your preliminary results much
for (i in seq(length(rawSubjectsData))) {
videoDataList[[i]] %<>%
dplyr::select(Leaves, Enters) %>%
apply(MARGIN=2, FUN=na.omit) %>% # removing NA observations
as.data.frame %>%
mutate(Enters = Enters - 600, Leaves = Leaves - 600) %>% # subtracting 600 seconds for warm-up
filter(Enters > 0 & Leaves > 0) # removing observations/tasks, which contained any of the warm-up period
}
Given the nature of the experimental study, the subjects performed varying tasks that may affect our gait analysis; we have only examined data where the subjects are outside of the main laboratory room since this captures the material handling component of the examined task. An immediate consequence of using only “out-of-the-room” data is that the time differences between the retained segmented gait cycles is no longer almost constant since relatively large time gaps exist, which correspond to when the participant is in the room. For this reason, we have subgrouped the data. Consecutive “out-of-the-room” gait cycles are placed in one subgroup. Statistically speaking, our subgrouping strategy follows the rational subgrouping mantra in statistical surveillance. The purpose of rational subgrouping is two-fold:
For more details on rational subgrouping, we refer the reader to Montgomery (2019) for more details. The code chunk below provides our approach to creating the rational subgroups for both the “profile” and “gait features” data sets.
numRows <- lengths(subjectProfilesList)/5
confidence <- videoDataList %>% lapply(function(x) {mean(x[,2] - x[,1])}) %>%
do.call(what=c) %>% `*` (0.05) # 5% confidence for each subgroup
subgroupList <- list()
for (id in seq(length(rawSubjectsData))) {
aM <- subjectProfilesList[[id]] %>% do.call(what=c) %>%
.[(numRows[[id]]+1):(2*numRows[id])] %>%
lapply(as.vector)
expT <- subjectProfilesList[[id]] %>% do.call(what=c) %>%
.[(4*numRows[[id]]+1):(5*numRows[id])] %>%
lapply(function(x) {x[1,1]}) %>% do.call(what=c)
###########################################################
# Match times of walking cycles, extracted from videos, against the experimental times
leaveIdx <- match.closest(videoDataList[[id]]$Leaves + confidence[id], expT)
enterIdx <- match.closest(videoDataList[[id]]$Enters - confidence[id], expT)
rmIdx <- c(3000, which(enterIdx - leaveIdx < 5))
leaveIdx <- leaveIdx[-rmIdx]
enterIdx <- enterIdx[-rmIdx]
########################################################## Create subgroups
aMList <- list()
LHDList <- list()
for (i in 1:length(leaveIdx)) {
aMList[[i]] <- aM[leaveIdx[i]:enterIdx[i]]
LHDList[[i]] <- subjectFeaturesList[[id]][leaveIdx[i]:enterIdx[i],]
}
######################################### Remove subgroups with length < 20
LHDList %<>% .[-{which(lengths(aMList) < 20) %>% c(3000)}]
aMList %<>% .[-{which(lengths(aMList) < 20) %>% c(3000)}]
######################################### Remove observations with profiles< 2nd percentile profile length; This is done to cut all profiles at this length (equal number of features)
cutLen <- aMList %>% do.call(what=c) %>% lengths() %>% quantile(probs=0.02)
for (i in 1:length(aMList)) {
rmIdx <- which(lengths(aMList[[i]]) < cutLen) %>% c(3000)
aMList[[i]] <- aMList[[i]][-rmIdx]
LHDList[[i]] <- LHDList[[i]][-rmIdx,]
}
######################################## Create a matrix with cutLen number of columns in each profile subgroup
aMList %<>% lapply(function(x) {do.call(rbind, x)}) %>%
lapply(function(x) {x[,1:cutLen]})
############################################ Lengths of subgroups 0&1: 10&20
subgroupList[[id]] <- list(Acc = list(X0=aMList[1:10], X1=aMList[11:30], X2=aMList[31:length(aMList)]),
LHD = list(X0=LHDList[1:10], X1=LHDList[11:30], X2=LHDList[31:length(LHDList)]))
}
names(subgroupList) <- c(paste0("subject", c(1:12, 14:15)))
Initialize.
T2Profiles <- list()
for (id in seq(length(subgroupList))) {
X0 <- subgroupList[[id]]$Acc$X0
############################# Phase0
############################# Remove single outliers in number of passes
pass <- 1
for (i in 1:pass) {
UCL0 <- BootstrapUCL0(x0=X0, alpha=0.05, iterations=200) # Find in Functions
######################### Find T2 for original X0 and remove outliers from subgroups
X0Pooled <- do.call(rbind, X0)
T2 <- X0Pooled %>% apply(FUN=mean, MARGIN=2) %>%
rep.row(n=nrow(X0Pooled)) %>% -(X0Pooled) %>%
apply(MARGIN=1, function(x) {x %*% ginv(cov(X0Pooled)) %*% x})
tmp <- mapply(function(x,y) {rep(y, x)}, x=sapply(X0, nrow), y=1:10) %>% do.call(what=c)
OCIdx <- T2 %>% split(f=tmp) %>% lapply(function(x) {c(which(x>UCL0), 3000)})
X0 %<>% mapply(FUN=function(x, y) {x[-y,]}, y=OCIdx) # Updated X0 after phase0
}
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- colMeans(X0Pooled)
############################# Within batch covariance; computed from phase0 after outlier removal
S_w <- X0 %>% lapply(colMeans) %>%
mapply(FUN=rep.row, sapply(X0, nrow)) %>% # Mean of each subgroup
mapply(FUN=function(x,y) {x-y}, y=X0, SIMPLIFY=F) %>%
lapply(function(x) {t(x) %*% x}) %>%
Reduce(f=`+`) %>%
`/` (sum(lengths(X0)/ncol(X0[[1]])) - length(X0))
############################# Phase1&2
############################# Read phase 1&2 subgroups
X1 <- subgroupList[[id]]$Acc$X1
X2 <- subgroupList[[id]]$Acc$X2
######################### Compute T2 in phase1
X1Mean <- sapply(X1, colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(X1)) %>% -(X1Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
######################### Bootstrap UCL using phase 1 data
sampled <- sample(T2, 1000*length(T2), replace=T) %>% matrix(ncol=length(T2), byrow=T)
UCL <- apply(sampled, 1, quantile, probs=0.99) %>% mean()
######################### T2 for phase2 data
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(c(X1, X2))) %>% -(X12Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
T2Profiles[[id]] <- list(T2=T2, UCL=UCL)
}
names(T2Profiles) <- paste0("subject", c(1:12, 14:15))
Initialize.
for (id in seq(length(subgroupList))) {
T2 <- T2Profiles[[id]]$T2
UCL <- T2Profiles[[id]]$UCL
color <- ifelse(T2>UCL, "OC", "IC")
T2.df <- data.frame(subgroup=1:length(T2), T2=T2, col=color,
shape=c(rep("Phase 1", 20), rep("Phase 2", length(T2) - 20)))
T2Plot <- ggplot(data=T2.df) +
geom_point(aes(x=subgroup, y=T2, shape=shape, col=col), size=2.5) +
geom_line(aes(x=subgroup, y=T2)) +
geom_hline(yintercept=UCL) +
annotate("text", 2, UCL, vjust=-1, label="UCL", size=5) +
scale_shape_manual(values = c(6, 16)) +
scale_color_manual(values = c("black", "red")) +
ggtitle(paste0("Subject ", id)) +
theme_bw() +
theme(plot.title=element_text(size = 20, hjust=0.5, margin = margin(t=10, b=-30)),
axis.text=element_text(size=15),
axis.title=element_text(size=20),
legend.title=element_blank(), # This removes the legend title
legend.position=c(.1, .68),
legend.text=element_text(size=10),
legend.key = element_rect(fill="grey"))
cat('###',paste0("Subject", id), "{-}",'\n')
print(T2Plot)
cat('\n \n')
}
Initialize.
# Reading the RPE values from gait_hotelling GitHub Repo:
RPE <- read.csv("https://raw.githubusercontent.com/fmegahed/gait_hotelling/master/Data/csvFiles/RPE/RPE.csv") %>% .[,-14]
SubgroupTimeList <- subgroupList %>% lapply(function(x) {x$LHD[c(2,3)]}) %>%
lapply(function(x) {do.call(c, x)}) %>%
lapply(function(x) {sapply(x, function(a) {median(a[,4])})})
OCIdx <- lapply(T2Profiles, function(x) {which(x$T2>x$UCL)})
for (id in seq(length(SubgroupTimeList))) {
OCTime <- SubgroupTimeList[[id]] %>% .[OCIdx[[id]]]
RPETime <- list(which(RPE[,id+1]>=13) %>% min(), which(RPE[,id+1]>=15) %>% min()) %>%
sapply(function(x) {RPE[x,1]})
cat('###',paste0("Subject", id), "{-}",'\n')
plot(OCTime, rep(1, length(OCTime)), pch=16, xlim=c(0,11000), xlab="Experiment Time",
ylab="", main=paste0("Subject", id))
abline(v=RPETime, col=c("blue", "red"), lwd=2)
legend("topleft", legend=c(expression(RPE >= 13), expression(RPE >= 15)), col=c("blue", "red"),
lty=1, lwd=2)
cat('\n \n')
}
Initialize.
T2Features <- list()
SWList <- list()
X0Updated <- list()
for (id in seq(length(subgroupList))) {
X0 <- subgroupList[[id]]$LHD$X0 %>% lapply(function(x) {x=x[,1:3]}) # Last column contains experiment time
############################# Phase0
############################# Remove single outliers in number of passes
pass <- 1
for (i in 1:pass) {
UCL0 <- BootstrapUCL0(x0=X0, alpha=0.05, iterations=200) # Find in Functions
######################### Find T2 for original X0 and remove outliers from subgroups
X0Pooled <- do.call(rbind, X0)
T2 <- rep.row(colMeans(X0Pooled), n=nrow(X0Pooled)) %>% -(X0Pooled) %>%
apply(MARGIN=1, function(x) {x %*% ginv(cov(X0Pooled)) %*% x})
tmp <- mapply(function(x,y) {rep(y, x)}, x=sapply(X0, nrow), y=1:10) %>% do.call(what=c)
OCIdx <- T2 %>% split(f=tmp) %>% lapply(function(x) {c(which(x>UCL0), 3000)})
X0 %<>% mapply(FUN=function(x, y) {x[-y,]}, y=OCIdx) # Updated X0 after phase0
}
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- colMeans(X0Pooled)
############################# Within batch covariance; computed from phase0 after outlier removal
S_w <- X0 %>% lapply(colMeans) %>%
mapply(FUN=rep.row, sapply(X0, nrow)) %>% # Mean of each subgroup
mapply(FUN=function(x,y) {x-y}, y=X0, SIMPLIFY=F) %>%
lapply(function(x) {t(x) %*% x}) %>%
Reduce(f=`+`) %>%
`/` (sum(lengths(X0)/ncol(X0[[1]])) - length(X0))
############################# Phase1&2
############################# Read phase 1&2 subgroups
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
######################### Compute T2 in phase1
X1Mean <- sapply(X1, colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(X1)) %>% -(X1Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
######################### Bootstrap UCL using phase 1 data
sampled <- sample(T2, 1000*length(T2), replace=T) %>% matrix(ncol=length(T2), byrow=T)
UCL <- apply(sampled, 1, quantile, probs=0.99) %>% mean()
######################### T2 for phase2 data
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(c(X1, X2))) %>% -(X12Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
T2Features[[id]] <- list(T2=T2, UCL=UCL)
X0Updated[[id]] <- list(X0=X0)
SWList[[id]] <- list(S_w=S_w)
}
names(T2Features) <- paste0("subject", c(1:12, 14:15))
names(X0Updated) <- paste0("subject", c(1:12, 14:15))
names(SWList) <- paste0("subject", c(1:12, 14:15))
Initialize.
for (id in seq(length(subgroupList))) {
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
color <- ifelse(T2>UCL, "OC", "IC")
T2.df <- data.frame(subgroup=1:length(T2), T2=T2, col=color,
shape=c(rep("Phase 1", 20), rep("Phase 2", length(T2) - 20)))
T2Plot <- ggplot(data=T2.df) +
geom_point(aes(x=subgroup, y=T2, shape=shape, col=col), size=2.5) +
geom_line(aes(x=subgroup, y=T2)) +
geom_hline(yintercept=UCL) +
annotate("text", 2, UCL, vjust=-1, label="UCL", size=5) +
scale_shape_manual(values = c(6, 16)) +
scale_color_manual(values = c("black", "red")) +
ggtitle(paste0("Subject ", id)) +
theme_bw() +
theme(plot.title=element_text(size = 20, hjust=0.5, margin = margin(t=10, b=-30)),
axis.text=element_text(size=15),
axis.title=element_text(size=20),
legend.title=element_blank(), # This removes the legend title
legend.position=c(.1, .68),
legend.text=element_text(size=10),
legend.key = element_rect(fill="grey"))
cat('###',paste0("Subject", id), "{-}",'\n')
print(T2Plot)
cat('\n \n')
}
Initialize.
# Reading the RPE values from gait_hotelling GitHub Repo:
RPE <- read.csv("https://raw.githubusercontent.com/fmegahed/gait_hotelling/master/Data/csvFiles/RPE/RPE.csv") %>% .[,-14]
SubgroupTimeList <- subgroupList %>% lapply(function(x) {x$LHD[c(2,3)]}) %>%
lapply(function(x) {do.call(c, x)}) %>%
lapply(function(x) {sapply(x, function(a) {median(a[,4])})})
OCIdx <- lapply(T2Features, function(x) {which(x$T2>x$UCL)})
for (id in seq(length(SubgroupTimeList))) {
OCTime <- SubgroupTimeList[[id]] %>% .[OCIdx[[id]]]
RPETime <- list(which(RPE[,id+1]>=13) %>% min(), which(RPE[,id+1]>=15) %>% min()) %>%
sapply(function(x) {RPE[x,1]})
cat('###',paste0("Subject", id), "{-}",'\n')
plot(OCTime, rep(1, length(OCTime)), pch=16, xlim=c(0,11000), xlab="Experiment Time",
ylab="", main=paste0("Subject", id))
abline(v=RPETime, col=c("blue", "red"), lwd=2)
legend("topleft", legend=c(expression(RPE >= 13), expression(RPE >= 15)),
col=c("blue", "red"), lty=1, lwd=2)
cat('\n \n')
}
Initialize.
qList <- list()
for (id in seq(length(subgroupList))) {
############### Read inputs
X0 <- X0Updated[[id]]$X0
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
X0PooledMean <- do.call(rbind, X0) %>% colMeans()
X1Mean <- sapply(X1, colMeans) %>% t()
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
S_w <- SWList[[id]]$S_w
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
t3 <- which(T2 > UCL)
p <- 3 # number of features
alpha <- 0.01
############## T2 Decomposition and UCL's
q <- init.decomp(p=p) # initialize the decomposition matrix
for(i in 1:nrow(q)) {# Control limits for each component of the decomposition
b <- subset(q[i,3:ncol(q)], q[i,3:ncol(q)] > 0)
T2 <- rep.row(X0PooledMean[b], length(X1)) %>% -(X1Mean[,b]) %>%
apply(MARGIN=1, function(x) x %*% solve(S_w[b,b]) %*% x)
sampled <- matrix(sample(T2, 1000 * length(T2), replace=T), ncol=20, byrow=T)
q[i,2] <- apply(sampled, 1, quantile, probs=1-alpha) %>% mean()
}
############## Decompose T2 for every out of control data point
tmpList <- list()
for(ii in 1:length(t3)){
for(i in 1:nrow(q)) {
b <- subset(q[i,3:ncol(q)], q[i,3:ncol(q)] > 0)
q[i,1] <- t(X0PooledMean[b]-X12Mean[t3[ii],][b]) %*%
solve(S_w[b,b]) %*%
(X0PooledMean[b]-X12Mean[t3[ii],][b])
}
############ Store decomposed elements
tmpList[[ii]] <- q
}
names(tmpList) <- paste0("subgroup ", t3)
qList[[id]] <- tmpList
}
names(qList) <- paste0("subject", c(1:12, 14:15))
Initialize.
################ Select subject number and features
id <- 1 # Select the subject number
f <- c(2,3) # Select 2 features by interpreting the T2 decomposition
################ Read inputs
X0 <- X0Updated[[id]]$X0
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- do.call(rbind, X0) %>% colMeans()
X1Pooled <- do.call(rbind, X1)
S_w <- SWList[[id]]$S_w
q <- qList[[id]][[1]] # Read one of the q's as its 2nd column (q[,2]) are the control limits and thus the same across all subgroups in one subject
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
################ Univariate CL's and 2D plot limits
XLim <- c(X0PooledMean[f[1]]-sqrt(S_w[f[1],f[1]]*q[f[1],2]),
X0PooledMean[f[1]]+sqrt(S_w[f[1],f[1]]*q[f[1],2]))
YLim <- c(X0PooledMean[f[2]]-sqrt(S_w[f[2],f[2]]*q[f[2],2]),
X0PooledMean[f[2]]+sqrt(S_w[f[2],f[2]]*q[f[2],2]))
xlim <- c(min(XLim) - 1.2*diff(XLim), max(XLim) + 1.2*diff(XLim))
ylim <- c(min(YLim) - 1.2*diff(YLim), max(YLim) + 1.2*diff(YLim))
################ Create grid for the ellipse contour plot
ngrid <- 100
grid <- cbind(seq(xlim[1], xlim[2], length = ngrid),
seq(ylim[1], ylim[2], length = ngrid))
################ T2 on the grid points
T2_grid <- apply(expand.grid(grid[,1], grid[,2]), 1,
function(x) {(x- X0PooledMean[f]) %*% solve(S_w[f,f]) %*% (x- X0PooledMean[f])})
T2_grid <- matrix(T2_grid, ngrid, ngrid)
######################################################## Plot
############## Plot the contour
# old.par <- par()
par(mar=c(5, 5, 4, 2) + 0.1)
plot(1 ,1, type = "o", xlim = xlim, ylim = ylim, main="Ellipse Chart", cex.main=1.5,
xlab=paste0(colnames(X0[[1]])[f[1]], " (m)"),
ylab=paste0(ylab=colnames(X0[[1]])[f[2]], " (s)"),
cex.axis=1.5, cex.lab=1.5)
legend("topleft", legend=c("phase 0 subgroup means", "phase 1 subgroup means",
"X control limits", "Y control limits"),
col=c("green", "blue", "black", "brown"), pch=c(3,3,NA,NA), lty=c(NA,NA,2,2),
cex=1, pt.cex=2, lwd=2, bg="grey")
################
idx <- apply(q[,3:5], MARGIN=1, function(r) all(c(f,0) %in% r)) %>% which()
contour(grid[,1], grid[,2], T2_grid, levels=q[idx,2], drawlabels=T, add=T)
################ Plot phase 0 data
center0 <- X0 %>% lapply(function(x) {x[,f]}) %>% sapply(colMeans) %>% t()
points(center0[,1], center0[,2], pch=3, cex=2, col="green", lwd=2)
# points(X0Pooled[,f], pch=1, cex=0.5, col="blue")
################ Plot phase 1 data
center1 <- X1 %>% lapply(function(x) {x[,f]}) %>% sapply(colMeans) %>% t()
points(center1[,1], center1[,2], pch=3, cex=2, col="blue", lwd=2)
# points(X1Pooled[,f], pch=1, cex=0.5, col="blue")
################ plot phase 1&2 data by subgroup number
center12 <- c(X1,X2) %>% lapply(colMeans) %>% sapply(function(x) {x[f]}) %>% t()
labels <- T2 %>% t() %>% as.data.frame() %>% names() %>% gsub(pattern='V', replacement='')
text(center12, labels=labels, cex = 1, col=ifelse(T2>UCL, "red", "black"))
################ Plot univariate CL's
abline(v=XLim, lty="dashed", col="black", lwd=2)
abline(h=YLim, lty="dashed", col="brown", lwd=2)
Baghdadi, A., Cavuoto, L. A., Jones-Farmer, A., Rigdon, S. E., Esfahani, E. T., & Megahed, F. M. (2019). Monitoring worker fatigue using wearable devices: A case study to detect changes in gait parameters. Journal of Quality Technology, Online First. https://doi.org/10.1080/00224065.2019.1640097
Maman, Z. S., Chen, Y.-J., Baghdadi, A., Lombardo, S., Cavuoto, L. A., & Megahed, F. M. (2020). A data analytic framework for physical fatigue management using wearable sensors. Expert Systems with Applications, 113405.
Maman, Z. S., Yazdi, M. A. A., Cavuoto, L. A., & Megahed, F. M. (2017). A data-driven approach to modeling physical fatigue in the workplace using wearable sensors. Applied Ergonomics, 65, 515–529.
Montgomery, D. C. (2019). Introduction to statistical quality control. Wiley.
Email: saebraga@buffalo.edu | Website: LinkedIn↩︎
Email: jiyeonk@buffalo.edu | Phone: +1-716-645-6063 | Website: University at Buffalo Official↩︎
Email: loracavu@buffalo.edu | Phone: +1-716-645-4696 | Website: University at Buffalo Official↩︎
Email: fmegahed@miamioh.edu | Phone: +1-513-529-4185 | Website: Miami University Official↩︎
Email: farmerl2@miamioh.edu | Phone: +1-513-529-4823 | Website: Miami University Official↩︎